# Direct outputs eligible for public release

# GLOBAL SETTINGS --------------------------------------------------------------

options(
    scipen = 999,
    digits = 16,
    max.print = .Machine$integer.max,
    show.error.locations = TRUE,
    warn = 1
)

RNGkind("L'Ecuyer-CMRG")
seed <- 818675309L
set.seed(seed) # setting main seed

# PACKAGES ---------------------------------------------------------------------
library(data.table)
library(ggplot2)
library(scales)

# PACKAGE SETTINGS -------------------------------------------------------------

# data.table
setDTthreads(threads = 1L)
options(datatable.print.class = TRUE, datatable.print.keys = TRUE)
# so that printing the data.table also shows the variable type on top

# BEGIN FILE -------------------------------------------------------------------

# READ IN OUR DATA -------------------------------------------------------------

moved_tract_dt <-
    readRDS(
        "~/estimation-output/wealth_effect_estimates_moved_tract.rds"
    )
setnames(moved_tract_dt, "ref_event_time", "event_time")
moved_tract_dt <-
    moved_tract_dt[
        model == "ratio" &
            event_time == 1 &
            income_quartile %in% c(1, 4, 5)
    ]
moved_tract_dt[income_quartile == 5, header_group := "Main"]
moved_tract_dt[income_quartile == 5, subgroup := "Full Sample"]
moved_tract_dt[income_quartile < 5, header_group := "Income"]
moved_tract_dt[income_quartile == 1, subgroup := "Quartile 1"]
moved_tract_dt[income_quartile == 4, subgroup := "Quartile 4"]
moved_tract_dt[, c("model", "income_quartile") := NULL]

moved_tract_subgroups_dt <-
    readRDS(
        "~/estimation-output/wealth_effect_estimates_moved_tract_subgroups.rds"
    )
setnames(moved_tract_subgroups_dt, "ref_event_time", "event_time")

moving_combined_dt <-
    rbindlist(
        list(
            moved_tract_dt,
            moved_tract_subgroups_dt
        ),
        use.names = TRUE
    )

# Rename some of the groups / subgroups, just for aesthetics
moving_combined_dt[subgroup == "Full Sample", subgroup := "Full\nSample"]
moving_combined_dt[subgroup == "Quartile 1", subgroup := "Quartile\n1"]
moving_combined_dt[subgroup == "Quartile 4", subgroup := "Quartile\n4"]
moving_combined_dt[
    header_group == "Presence of Young Children",
    header_group := "Presence of\nYoung Children"
]
moving_combined_dt[
    header_group == "Age of Winner",
    header_group := "Age of\nWinner"
]
moving_combined_dt[
    subgroup == "Younger Winner",
    subgroup := "Younger\nWinner"
]
moving_combined_dt[
    subgroup == "Older Winner",
    subgroup := "Older\nWinner"
]
moving_combined_dt[
    header_group == "Strong Work Attachment",
    header_group := "Strong Work\nAttachment"
]
moving_combined_dt[
    header_group == "Homeownership",
    header_group := "Home-\nownership"
]

moving_combined_dt[
    ,
    header_group :=
        factor(
            header_group,
            levels =
                c(
                    "Main",
                    "Income",
                    "Strong Work\nAttachment",
                    "Presence of\nYoung Children",
                    "Age of\nWinner",
                    "Home-\nownership"
                )
        )
]
moving_combined_dt[
    ,
    subgroup :=
        factor(
            subgroup,
            levels =
                rev(
                    c(
                        "Full\nSample",
                        "Quartile\n1",
                        "Quartile\n4",
                        "No",
                        "Yes",
                        "Younger\nWinner",
                        "Older\nWinner"
                    )
                )
        )
]


figure_5_3 <-
    ggplot(
        data = moving_combined_dt,
        aes(
            x = subgroup,
            y = estimate * 100000 * 100
        )
    ) +
    geom_bar(stat = "identity", position = position_dodge(), fill = "black") +
    geom_errorbar(aes(
        ymin = ((estimate - (1.64 * cluster_se)) * 100000 * 100),
        ymax = ((estimate + (1.64 * cluster_se)) * 100000 * 100)
    ),
    position = position_dodge(0.9),
    show.legend = FALSE,
    width = 0.2,
    color = c("#999999")
    ) +
    theme_bw(base_size = 10) +
    scale_y_continuous(
        breaks = (seq.int(from = 0, to = 1000, by = 50) / 100),
        labels = scales::number_format(
            accuracy = 0.1,
            big.mark = ","
        )
    ) +
    labs(x = NULL, y = "Effect per $100,000 of Lottery Winnings (pp)") +
    coord_flip(ylim = c(0, 7)) +
    facet_grid(
        header_group ~ .,
        space = "free",
        scales = "free",
        switch = "y"
    ) +
    theme(
        panel.grid.minor = element_blank(), # remove actual gridlines
        panel.grid.major = element_blank(),
        legend.position = "bottom",
        legend.title = element_blank(),
        legend.key.size = unit(1.75, "lines"),
        strip.placement = "outside",
        panel.spacing.y = unit(0, "cm"),
        strip.text.y = element_text(angle = 270, size = 6.5),
        panel.border = element_blank(),
        axis.line = element_line(),
        axis.text.y = element_text(hjust = 0.5)
    )

ggsave(
    plot = figure_5_3,
    filename = "~/paper/figures/figure-5.3.png",
    width = 6,
    height = 4.5,
    dpi = 600
)

ggsave(
    plot = figure_5_3,
    filename = "~/paper/figures/figure-5.3.tif",
    width = 6,
    height = 4.5,
    device = "tiff",
    dpi = 600
)

# Housekeeping
moved_tract_dt <- NULL
moved_tract_subgroups_dt <- NULL
moving_combined_dt <- NULL
figure_5_3 <- NULL
rm(
    moved_tract_dt,
    moved_tract_subgroups_dt,
    moving_combined_dt,
    figure_5_3
)